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Figure 9 : Locations of Monitoring sites for Ambient Sulphate Concentrations 
(Marked by Numbers) and Wet Deposition (Marked by Letters) 
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Figure 8 : Model Computed Annual Accumulation, wet Deposition of Total Sulphur 
in kg ha yr~ , Including Background of 2 kg ha"~ 1 yr~ 1 
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Figure 7 : Model Computed Annual Average S0 4 Concentration for 1 979 
In jig m 3 Including Background of 2/Lig m~ 3 
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Figure 6 : Model Computed Annual Average SO Concentration for 1979 in/jg m 3 
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Figure 5 : Model Computed Monthly Accumulation, wet Deposition of 

Total Sulphur for July, 1979 in Units of 10/xmole m-2Month~ 1 
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Figure 4 : Model Computed Monthly Average S0 4 Concentration for 
July, 1979 inug m" 3 , Including Background of 2yg m" 3 
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Figure 3 : Model Computed Monthly Average S0 2 Concentration 
for Jury, 1979 in fxg m-3 
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Figure 1 : Model Emission Inventory for 1979. Size of Circle Indicates Source Strength 



TABLE 5 



Wet Deposition Total Sulphur 



Year 1979 



«fj 



Site 



Name 



Model 


Observed 


Value 


Value 




(prorated) 


kg ha" 1 y" 1 


kg ha" 1 y" 1 



Whiteface Mt. 


A 


5.7 


Ithica 


B 


7.6 


Penn. State 


c 


9.9 


Oxford Oh. 


G 


5.9 



MAP3S Network 
1979 

7.1 

8.7 

9.7 
10.4 



Ela - Kenora 


H 


2.9 


Long Point 


I 


13.3 


Chalk River 


J 


6.7 



APN Network 

1979 

3.2 
13.3 
11.1 



Atikokan 


K 


3.0 


Kingston 


L 


6.4 


Peterborough 


M 


10.2 


Simcoe 


N 


13.3 



CANSAP Network 
1978 

3.5 

9.1 

8.5 
13.9 



background of 2 kg ha 1 yr~ included. 



GE/mb/AQM9-13 



TABLE 4 



SO, Concentration 
Year 1979 



te- 



station 



Computed 



** 



Measured 



Corrected 
(-5 ug m" 3 ) 



Ratio 





ug m- 


ug m 


ug 


m 


Computed 










Measured 






MOE (1979) 








1. Sudbury 


9.5 


10.1 


5. 


.2 


1.8 


3. Toronto 


5.8 


13.5 


8.5 


0.7 


4. Hamilton 


6.8 


15.0 
SURE* (1978) 


10, 


.0 


0.7 


5. Detroit 


9.1 


7.0 






1.3 


6. Indiana 


6.8 


7.5 






0.9 


7. Indiana 


8.0 


6.9 






1.2 


8. Indiana 


8.0 


7.8 






1.0 


9. Ohio 


6.7 


5.7 






1.2 


10. Ohio 


7.3 


8.8 






0.8 


11. Ohio 


9.1 


7.5 






1.2 


12 Penn. 


6.7 


6.7 






1.0 


13 Penn. 


9.6 


8.0 






1.2 


14. Penn. 


5.3 


7.4 






0.7 


15. Penn. 


6.6 


7.4 






0.9 


16. Illinois 


9.3 


5.6 






1.7 


17. Illinois 


6.0 


3.3 






1.8 


18. Illinois 


6.6 


6.0 






1.1 


19. New York 


5.7 


5.3 






1.1 


20 New York 


5.5 


6.1 






0.9 


21 New York 


4.4 


5.4 






0.8 



*.,. 



Nieman analysis Aug. 1977 to June 1978 



** _3 

Includes background of 2 ug m plus background from 21 additional sources 

as calculated by statistical model. 



TABLE 3 
Wet Deposit on of Total Sulphur 

July 1979 



Site 



Name 



# 



Model 




Observed 


Value 




Value 


umole 


m 


-2 -1 
mo 











MAP3S Network 


Whiteface Mt. 


N.Y. 


A 


701 


2010 


Ithica 


N.Y. 


B 


1960 


5340 


Penn. State 


Penn. 


C 


4209 


3600 


Charl. 


Va. 


D 


2142 


1710 


Urbana 


111. 


E 


3970 


4140 


Lewes 


Del. 


F 


2360 


2520 


Oxford 


Oh. 


G 


4926 


7,700 
APN Network 


Ela-Kenora 


Ont. 


H 


32 


300 


Long Point 


Ont. 


I 


4955 


4,000 


Chalk River 


Ont. 


J 


1049 


4,000 



TABLE 2 



SO. CONCENTRATION 
4 






Site 



Name 



JULY 


1979 




Model* 




Observed 


Value 




Value 


-3 

ug m 


Reading 

-3 
ug m 


Corrected 
(-5 S g~ 3 ) 
ug m 



Sudbury 
Parry Sound 
Toronto 
Hamilton 



Detroit 
Indiana 



Ohio 



Penn. 



Illinois 



New York 







MOE Network 






July 1979 


1 


7.8 


12 


2 


5.2 


11 


3 


6.6 


16 


4 


8.4 


16 

SURE Network 
July 1978 


5 


9.3 


8 


6 


5.1 


11 


7 


8.7 


12 


8 


8.3 


13 


9 


6.8 


14 


10 


7.0 


14 


11 


15.4 


17 


12 


8.0 


13 


13 


15.5 


16 


14 


7.1 


12 


15 


9.5 


11 


16 


13.9 


10 


17 


11.1 


8 


18 


10.0 


15 


19 


6.6 


15 


20 


6.4 


10 


21 


5.3 


10 


22 


3.8 


7 



7 

6 

11 

11 



* -3 

Background of 2 ug m included. 



TABLE 1 
Values of Model Parameters 



July 1979 Year 1979 

SO2 to SO4 transformation rate 2% hr _1 1% hr -1 



v d 


SO2 dry deposition velocity 


1 cm s" 1 


vd 


SO4 dry deposition velocity 


0.1 cm s _1 


H 


mixing height 


1000 m 


d 


initial deposition of SO2 





6 


local oxidation SO2 to SO4 


10% 


At 


time step 


6 hour 


k w 


SO2 wet depositon rate 


(3 x 10-^)1 


w 


SO^= wet deposition rate 


(3 x 10-<0I 


1 


precipitation rate 


mm s"l 



5% 



■_, 
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* - 
6. Conclusion 

The Lagrangian model which has been described in this report 
constitutes a significant departure from the statistical model, in that 
the motion and the state of pollutant is calculated step by step during 
transport along trajectories. Over a long period of time these 
trajectories may be considered collectively as a statistical ensemble. 
The mean trajectory and standard deviation may then be regarded as 
the centre line of the path of the pollutant and the large scale 
horizontal diffusion. This is in contrast with the statistical model 
which assumed these values at the outset. Where the statistical 
model assumed single values for calculating mass decay and 
conversion, in the Lagrangian model an attempt is being made to 
simulate temporal and spatial changes in these factors. The 
Lagrangian model consumes much more computing time and hence is 
a great deal more costly to run. Because of the interactions among 
the various components, the results of changes in the model cannot 
be anticipated with certainty. So sensitivity tests are necessary 
before incorporating these changes. 

Such influences as occurrence of rain, small variations in the 
winds and the chemical changes in the pollutant cannot be defined 
precisely at each step. This means that the trajectories, 
concentrations and depositions calculated by the model must be 
averaged over a period of time to be considered valid. 

The basic form and mechanisms of calculation of the model are 
now satisfactorily developed. However, considerable testing and 
adjustments are necessary before model predictions can be regarded 
with confidence. 



in the lowest 50 m. Near the surface the ratio C /C = 0.6 when 
V<j=l cm s -1 . Graphs by Scriven and Fisher (1975) of C vs distance 
from the source, for values of V^ from 0.1 to 5.0 cm s** show 
modified exponential decay of C with distance. An almost linear 
decrease of C with increasing values of V<j illustrates the 
compensating relationship between the two parameters. 

5 A Transformation of SO? Sulphur dioxide to Sulphate 

Photochemical effects on the rate of conversion can be 
represented by diurnal and seasonal changes in k-^ or by a dependence 
on cloud cover and the sun's angle. The conversion rate appears to 
increase in rain. Thus k t could be increased for wet sections of the 
trajectories. Transformation also has a dependence on the 
concentration of SO^~ and this is more difficult to calculate because 

it makes eq. (1) non-linear. It may be possible to introduce an 
iterative calculation for this process. This can be tried, obtaining a 
first value of concentration, then running the model repeatedly with 
k t adjusted to the previously obtained concentration each time. 






and boundary conditions 

_JLi£s 2H_ = atz=H (12a) 

di 

Ik- ^ 
_Jl£i £*__ =V d .C atz=0 (12b) 

For SO2 the SINK terms in equation (11) express deposition (dry 
and wet) and transformation to SOjp; for SO^ = the SINK terms are 
the depositions only. The point source for SO2 is located either at 
ground level or elevated a fraction of the boundary layer height. 
There are two source terms for SO^=, the point source and the source 
due to chemical conversion from SO2 and hence proportional to the 
transformation loss term of the SO2 equation. 

Shannon (1981) integrates equation (11) numerically, using 
several vertical levels. r . Fisher (1975) and Maul (1977) have solved 
equation (11) analytically for constant k z over the boundary layer, 
taken either as a whole or in parts. The near surface values of k z are 
of primary importance in the calculation of deposition. Above 50- 
100m, k z can be considered a constant for long range transport 
calculations. Either a numerical solution such as Shannon's or an 
analytical solution could be used in the model to simulate the vertical 
concentration distribution, with attention directed to the near- 
surface values. As previously noted, it is important to test the 
sensitivity of the model to the range of possible values of C before 
embarking on this type of complication in the model process. A 
simple adjustment in the value of V^ may be found to be sufficient. 

Graphs of variation of concentration with height as obtained 
from the analytical solution by Maul (1977) show that beyond 100 km 
the concentration is distributed almost evenly in the vertical, except 



':" J 



atmospheric stabilities. A survey of the literature by Sehmel (1980) 
shows that values of V<j measured by many different authors is within 
this range, with some extreme values ranging as low as 0.0k and as 
high as 7.0. Values for snow surfaces may have to be added to Sheih's 
tables which assume midsummer conditions. These have been 
measured by Whelpdale and Shaw (1974) and Barrie and Walmsley 
(1978). 

5.3.2 Mixed Layer Depth 

The height to which surface mixing is effective rises in the 
daytime when the heated surface warms the lower atmospheric layers 
and lowers at night when the surface cools. This effect can be 
included in the model by diurnal variation of H. A more explicit 
calculation of H can be based on values of stability calculated from 
cloud cover and the sun's angle. When a trajectory moves into a 
region of lower H value, some of the mass of the pollutant remains 
above the boundary layer. The model will then suffer loss of mass at 
the upper boundary due to this effect, unless the mass is recaptured 
when the value of H is higher. 



5.3.3 Surface Concentration, Vertical Dispersion 

When dry deposition is expressed in terms of surface 
concentration C , as opposed to the layer-average value C, then the 
time rate of change in concentration is expressed by the diffusion 
equation in the vertical, 



cX 



\ 



CJ KV_ c)£- + SOURCES + SINKS (11) 

air 
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5.3 Dry Deposition 

The calculation of dry depsition for both SO2 and SO^ = can be 
generalized for the purpose of this discussion as 

Drydep = k d .C = V d .c", (10) 

H 
where C is the layer averaged concentration. Both Vj and H 
currently are assigned fixed values. 

In reality V d and H both vary with thermal stability and the 
type of surface, e.g. vegetation type and water surface. Also the 
appropriate concentration is not C but the surface value C . 

It should be noted here that there is a tendency for the model 
to be insensitive to changes in dry deposition. This is because of the 
compensating relationship which exists between k d and C. If kj is 
made smaller, then downstream values of C are greater. Hence the 
product k d .C does not change as much as the terms taken 
individually. Consequently complicated calculations may not be 
justified in the search for representation of reality here. 

5.3.1 Deposition velocity 

Sheih et al (1979) have produced maps of variations is V d based 
on land use. These values could be incorporated into the model. The 
overall range of values for V d is from 0.1 to 1.6 cm s~*. These values 
were obtained from field measurements for varying surface types and 



(3) The problem of the "patchiness" of rain arises because of the 
variability of rain on time and space scales which are smaller than 
that of the model. Rain clouds are meso-scale (~10 km) phenomena 
which cannot be resolved on the 127 km grid. In addition, rain 
occurrence varies on a time scale of about 15 minutes, whereas the 
precipitation effect in the model is calculated at 6 hour intervals. 

Smith (1981) attempted to solve this problem using statistical 
methods. For example, for a given location, he determined the 
probability that there will be no rain in the next 15 minutes, when 
rain has occurred in the preceding 6 hours. In addition to this type of 
calculation, he separated the trajectories into wet regions and dry 
regions, computing statistics of the movement from wet to dry 
regions and vice versa. 

It would appear then that the use of a probability factor p r 
might be introduced into the wet deposition computation, making the 
general formula 

Wet dep = p r (WI) C (9) 

H 
where C is the concentration, i.e., q or s. Determination of p r would 
take into account, in some way, the considerations given above. 
Concerning the movement of trajectories from wet to dry regions, it 
may be useful to examine the relative motions of rain clouds and 
trajectories, by expressing the correlation between boundary layer 
and cloud level winds. Thus if a trajectory is in a wet region there 
may be some tendency for it to move with the rain area, thus 
increasing the value of p r . 



5.2 Wet Depostion 

Three types of modification are suggested here. These are (1) 
subdivison of a trajectory segment into two parts if it passes from 
one grid square to another during a time step; (2) simple 
parameterization of cloud chemistry, (3) incorporation of a 
representation of the "patchiness" of rain (Smith; 1981). 

(1) When a trajectory segment crosses from one grid square to 

another, the calculation of wet deposition can be broken into 

two parts using the appropriate k^ and k w values for the two 

grid squares. 

_ ' '' ■ •. 

(2). The parameter 1^ (or k w ) can be expressed (Barrie; 1981) as 



kw =WI 
H 



where W is a washout coefficient dependent on temperature and the 
pH of rainwater. To use this in model calculations, it would be 
feasible to designate typical seasonal temperatures and obtain the 
appropriate values of W from Barrie's graph of W vs pH for given 
temperatures. It would be more difficult, however, to define 
appropriate pH values in the model. A single representative value 
could be assumed, as a first approximation. 



5.1 Wind in the Boundary Layer 

Instead of varying the wind speed and turning angle diurnally to 
parameterize stability effects as is now being done, it is possible to 
use explicit values of stability for this purpose. The simplest method 
is to define stability categories based on values of insolation 
computed from cloud cover data and the angle of incidence of the 
sun's radiation. 

A more complex method of computing wind values directly 
using similarity principles was devised by Danard (Danard et al; 
1980). It is comparable in theory to the calculations of Peterson et al 
(1981, pp. 17-27) for the Windatlas of Denmark, but differs in 
practical application. Danard's technique is designed for modelling. 
It employs an iterative procedure for computing stability parameters 
using graphs derived by Yamada (1976) from measured values. Input 
information used by Danard are the wind and temperature at the top 
of the boundary layer plus, if possible, surface temperature. If this is 
not available a second iteration scheme is linked to the first. In it 
the surface temperature is computed from an equation of the energy 
balance at the earth's surface. Information on incident solar 
radiation and humidity are needed for this calculation. The two 
iteration procedures are mutually dependent. 

Some adaptation of Danard's method may be possible using 
surface geostrophic winds computed by the model and estimates of 
temperature profiles which are typical seasonally. 



in addition, are difficult to formulize in a simplified manner. Hence 
their effects on model results should be assessed before adopting 
then. 

After the model has been tested as described above, it can be 
applied to data now being accumulated on the Eclipse computer for 
1981 and later. 

Some efforts can be directed also to calculating atmospheric 
transport of nitrates. The source inventory will have a different 
form in that the nitrogen emissions are mostly from vehicles on 
highways. In addition, the atmospheric chemistry of nitrogen is 
probably less well simulated by linear interactions. However this can 
be attempted as a first approximation. 

Finally, it has been suggested that it may be possible to design 
an advanced statistical model based on the form of the present 
statistical model, using trajectory statistics from the Lagrangian 
model. A consideration of this idea will be given at a later date. 
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5. Future Modifications of the Model 

After the calculations for 1979 and 1978 have been completed 
there should be a review of the model results to assess the strengths 
and weaknesses in the model simulations. 

Since the trajectory calculations form the basic unit of the 
model, it would be well to examine these first. Comparison with 
examples of 850 mb. trajectories can provide insight into the relative 
merits of using surface and upper air data for these computations. 
Calculation of the deviation of the boundary layer wind from 
geostrophic values due to surface friction and stability effect merits 
further consideration. Possible modifications to these calculations is 
discussed in section 5.1. 

Having the collection of trajectories as a base, the mass 
balance calculations can be examined to determine how 
parameterization can be improved. Wet deposition is of primary 
concern, because of its influence on concentration calculations and 
because of the difficulties in parameterizing its effect. These 
problems arise from cloud physics and chemistry and because rain 
occurs unevenly on spatial and temporal scales which are not resolved 
by the model. Possible modifications to wet deposiiton calculations 
are outlined in section 5.2. 

Dry deposition and the changes in its parameters are discussed 
in section 5.3. Possible changes in the calculation of transformation 
of SO2 to SO^" are given in section 5A. 

Tests of sensitivity of the model to changes in the parameters 
should be performed before any new complexities are introduced. 
The parameters are frequently based on chemical and physical 
influence in the atmosphere which are not well understood and which 



Underprediction occurs at Oxford Oh. where the observed value is 
quite high. As in the July observations, this may result from a 
localized source. 



compared to measured values. The SURE values from U.S. monitors 
were analysed by Niemann (1981) for the period August 1977 to June 
1978. Although this is a different time period than the modelling 
period, the values from the Niemann analysis should give an 
approximation of the concentrations for 1979. MOE measured data 
had a correction factor of -5 pg m~3 applied as for the July data. In 
the MOE measured data, Toronto and Hamilton reported daily more 
than 70% of the year. Sudbury reports were daily for July and 
August, but only every six days for the rest of the year, with a 
sample size of 76 days for the year. Measured annual concentrations 
are lower than for the month of July. This might be expected since 
SO^ concentrations usually peak in summer. Computed annual 
concentrations are more evenly distributed from site to site than 
they were for July alone. This is due to the variations in the many 
trajectories computed for the full year. Good agreement is shown 
between computed and observed values. The range of the Ratio: 
Computed/Measured is from 0.7 to 1.8. Near the source areas, the 
concentration is overpredicted. This is evident in the Ratio column 
where #1 Sudbury has a value 1.8, #5 Detroit, 1.3, #16 Illinois, 1.6 
and #17 Illinois, 1.8. All other ratios are near 1 or less than 1, the 
lowest being at Toronto and Hamilton with values of 0.7. 
Wet Deposition 

The computed pattern of wet deposition of total sulphur in Fig. 
8 shows a broad maximum in the Ohio valley extending into 
southwestern Ontario, with a smaller maximum near Sudbury. 
Comparison with measured values in Table 5 shows good general 
agreement between observed and model values. At sites distant from 
the major sources areas, e.g., Whiteface, Ithica, Ela-Kenora and 
Chalk River, wet deposition is underpredicted as in the July case. 
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4.3 One year simulation 

The first tests of the model for a full year's data gave results 
which are discussed below. These must be regarded as preliminary 
results since there has been no tuning of the model to a full year of 
meteorological data. With the relatively large number of trajectories 
for a year, their variations provide horizontal dispersion. So no 
additional dispersion was applied such as the Gaussian distribution for 
the one month simulation. 

4.3.1 Model results 

The computed annual concentration of SO2 is shown in Fig. 6. 
Two maxima exceeding 40 ug m - ^ occur in the Ohio valley near the 
source regions and a maximum also occurs at Sudbury. The 5 ug m~3 
isopleth encloses much of the central area of the northeastern states, 
extending into Ontario. 

Fig. 7 shows the computed annual concentration of SO4. The 
background value of 2 ug m"^ has been added. This represents a 
generally accepted value of ambient sulphate from natural sources 
plus small distant anthropogenic sources not included in the emission 
inventory. In addition, concentrations and wet depositions resulting 
from the 21 sources not included in the model were computed by the 
statistical model and added. As may be expected, the SO4 

concentrations do not show the steep gradients that the SO2 isopleths 
exhibit. The smooth variation of SO4 values is in agreement with 
observations (see Table 5). Maxima in excess of 10 jig m~3 are 
computed for the source areas. The 6 ^<j m~3 isopleth encloses the 
central northeastern states and extends into southwestern 
Ontario, with a 6 ug m"^ maximum near Sudbury. 

In Table 4, computed annual sulphate concentrations are 



.:: > 



Wet Deposition 

Accumulated wet deposition of total sulphur for the month as 
computed by the model is plotted in Fig. 5. Maximum values 
exceeding 4,000 u mole m"2 occur in the states adjacent to the major 
source areas. The location of this area of maximum deposition 
appears to be in agreement with reality (MAP3S/RAINE; 1980). Wet 
deposition exceeding 1,000 u mole m~ 2 is computed for the wider 
area of all the central states of the northeastern U.S., all of southern 
Ontario, projecting into northern Ontario and Quebec. Table 3 shows 
computed wet deposition of sulphur compared with amounts observed 
at MAP3S monitoring sites for the month of July, 1979 
(MAP3S/RAINE: 1980) in the U.S. and the APN network (Barrie; 1981) 
in Canada. Where computed values differ signficantly from observed, 
they tend to be low, especially at sites which are far from the source 
areas such as Whiteface Mt. (ft A), Ithica ( #B) and Ela-Kenora ( # H). 
As in the case of the concentration the scarcity of trajectories 
through distant grid squares for a single month simulation is one 
contributing factor to these low values. Another is the omission of 
sources outside the grid, e.g. east of New York State and west of 
northern Ontario. Wet depositon is overpredicated at Penn State (# 
C), Charlottesville (# D) and Long Point (# 1^ near the source areas. 
Oxford ( #G) shows a high predicted wet deposition, but not as high as 
observed. However, the observed value is unusually high at this site. 
Possibly this is due to a very localized source, (*v50 km) which the 
model cannot duplicate. 



Table 2 compares SO^= concentration values computed by the 
model with measured values for specific sites monitored by the MOE 
network in Ontario and the SURE network in U.S. The locations on 
the grid array were chosen by eye and are judged to be within + ft 
grid distance of each site. Locations are shown in Fig 9. A constant 
background of 2 ug m -* has been added to computed concentrations. 
Recorded SO^ = concentrations from the MOE network were 
corrected for spurious high readings which result from conversion of 
SO2 to SO^ = on the Gelman AE glass filters used in the HIVOL 
monitors. The error is estimated by Lusis (1981) as +6 to +7 ug m~3 
and by Stevens (1981) as +3 to +4 ug m~3. Based on these findings a 
constant factor of 5 ug m"3 was substracted from all MOE monitored 
values, as shown. SURE measured data were taken from an analysis 
by Niemann (1980). These are for July 1978 and hence are not 
directly comparable to computed July 1979 values. But they should 
provide some indication of the appropriate concentrations. 

The model tends to overpredict SO4 concentrations near the 
sources, such as at //l Sudbury, #5 Detroit and #16 and #17 in Illinois. 
Away from the source areas, concentrations are underpredicted. This 
may be due in part to localized small sources, possibly an important 
factor for urban monitoring sites such as Toronto and Hamilton. One 
of the reasons for the low computed values in New York state is the 
exclusion from the model of sources just to the east in neighboring 
states. Also, if the trajectory is dry, i.e, no precipitaiton scavenging 
occurs, there is likely to be sulphur still in the air at the end of 48 
hours. This suggests that, if trajectories were to be extended to a 
longer period, e.g. 72 hours, the sulphur would be transported farther 
to the points remote from the sources. 






where Kq is the diffusivity coefficient. The value Kq =10^ m^s - * 
simulates sub-grid wind fluctuations on the meso-scale plus the 
effect of vertical wind shear in creating horizontal wind variations at 
different levels in the vertical, as discussed by Smith (1979). This is 
intermediate between values of Krj = .0^ and lO^m^s - * derived 
empirially for synoptic scale trajectory variations by Murgatroyd 
(1969) and Durst et al (1959), and 10^ m 2 s" 1 , employed by 
Bhumralkar et al (1980) for meso-scale diffusion. For Kq = 10* m 2 
s-* the puff is spread over a width of approximately 3 grid squares at 
t=6 hours and 9 grid squares at 48 hours. 

4.2.2 Model results 
Concentration 

In Figs. 3 and 4 which are plots of computed SO2 and SO^ = 
concentration, high values of SO2 and SO^= concentrations generally 
correlate with source areas. Due to the large value of k^ for SO2, 
the ambient concentration of SO^ = in the model is influenced 
strongly by its initial value. Thus this result is not surprising. Two 
maxima of 30 ug m~ 3 of SO2 occur in the Ohio River 

valley. Near Sudbury Ontario the model computes a maximum SO2 
concentration of 20 ug m -3 . The 5 ug m _ 3 line encloses an area 
encompassing most of the central northeastern states and the area of 
southern Ontario close to the lower Great Lakes. 

Concentration of SO^ = varies more smoothly than S02» because 
of the continuous conversion along the trajectories. The two SO^= 
concentration maxima in the Ohio River valley are 10 and 15 ug m - ^. 
A maximum in the Sudbury area exceeds 6 ug m~3. The 3 ^jg m~3 
line encompasses the northeastern states, southern Ontario and 
reaches northward to the area around Sudbury and into Quebec. 



model, which is designed for long term calculations (1-2 years) is 
applied to a shorter time period such as one month. These 
modifications include a change in the value of the initial conversion 
<j>and of the hourly conversion rate k t as discussed in section 3.1. In 
addition, the introduction of horizontal dispersion calculations is an 
important modification for a one month simulation. 

4.2.1 Horizontal Dispersion 

Because of the relatively small number of trajectories in a one 
month simulation (about 120 per source for 6-hour intervals), the 
mass of pollutant is not distributed sufficiently broadly by the 
variation in trajectory paths. As a result, concentrations show large 
discontinuities between neighboring grid squares. For the purpose of 
smoothing, horizontal dispersion has been applied, expanding each 
puff in time as it moves along the trajectory. The result is a set of 
'fuzzy' trajectories. The use of 'fuzzy' trajectories enables us to 
obtain better averages than would be obtained from individual 
trajectories. Ideally one should average over several years for the 
same month. The 'fuzzy' trajectory reflects approximately the 
variability of a single trajectory from year to year. Horizontal 
dispersion spreads the centre line concentration and deposition values 
over several grid squares. A normal or Gaussian distribution has been 
assumed for each puff such that the concentration C at a point x, y 
and time t is related to the centre line mass Q' by the equation 



C(x,y,t) = Q' (x,y,t) l*exp 

H 2-ir^<T y 



-(x-a 2 - (y-y)2 
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Assuming (T x = (Ty = (T» the puff radius is 2 07 
Assigning an inital value <T Q = 0.25 grid interval, 
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values derived from theoretical considerations. These are subject to 
continuous review and testing. Laboratory measurements do not 
always transfer well to the field. On the other hand field 
measurements are subject to many difficulties and are very expensive 
to carry out. Also, values derived under one set of meteorological 
conditions may be appropriate only to those conditions. Thus 
parameter values are not well known. Modelling makes use of results 
of research in representing physical and chemical interactions. It 
also contributes to an understanding of these processes by providing a 
different type of experimentation in which the numerous parameters 
may be tested. 

4.1.3 Nitrogen oxides 

Use of the model has been proposed to simulate the transport of 
the oxides of nitrogen. These pollutants undergo more complex 
chemical reactions in the atmosphere than sulphur oxides. Hence the 
use of linear equations in the model is even less applicable when 
applied to sulphur. It may be possible to apply it as a first 
approximation, with subsequent suitable refinements. 

4.2 One Month Simulation 

Average ambient SO2 and SO^= concentrations and wet and dry 
depositions have been computed for the month of July 1979. This is 
for the 62 emission sources within the grid. Calculations and results 
for the single month's data are described in the following sections. 
These are preliminary results. Considerable testing with variations in 
the model parameters will be necessary before the model can be 
considered valid. 

Modifications of model calculations are necessary when the 



scale with tests on seasonal, monthly and possibly episode (about 4 
days) scale. 

Planned acquisition of an improved meteorological data set for 

1979 and for 1980 plus use of the meteorological data being stored 
continuously on the Eclipse computer beginning in the latter half of 

1980 will permit simulations for the years 1979, 1980, 1981 and 
future years. 

Some adaptations of the model to the 1978 data are 
necessary. An increased grid area is possible because of improved 
geographical coverage in the meteorological data. Data is hourly 
from Canadian stations and 3-hourly from U.S. stations, as opposed to 
the 6-hourly data used in the calculations for 1979. For 1978 
precipitation data is also different. It is descriptive rather than 
quantitative. Precipitation is classified as to rain or snow, 
continuous or showery and is described as light, medium or heavy in 
these categories. 

Measured values of concentrations and depositions are 
necessary for comparison with model results in order to evaluate the 
model's performance. These are being assembled. They are obtained 
from many diverse measuring networks in Canada and the United 
States. These networks have developed independently and for 
different purposes. In consequence the types of data and their 
presentation are quite varied. It will take time to organize their data 
for comparison with the model. However, preliminary comparisons 
can be made using parts of the data. 

4.1.2 Model Parameters 

Testing and modification of mass balance parameters is a 
continuing process in model development. The parameters values are 
taken from the literature. They include results of measurements and 
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4. Present State of the Model and Future Projections 

Both major elements of the model, i.e. trajectory and mass 
balance calculations have been developed. Results of preliminary 
tests indicate that the mechanics of the model are working 
satisfactorily and that predicted concentrations and depositions are 
the right order of magnitude. Proposed work for the coming months 
is outlined in section 4.1. 

The simulation of one month's data for 3uly 1979 has been used 
as a base for testing the model parameters. Discussion of the July 
1979 simulation is presented in section 4.2. All the trajectories have 
been calculated for the full year 1979. They have been combined with 
the gridded precipitation data and stored in such a way that the 
precipitation data accompanies each trajectory segment. A 
preliminary simulation of the mass balance for the year has been 
completed. Results from this first test are presented in section 4.3. 

4.1 Proposed Working Schedule 

A graphical representation of a proposed working schedule is 
shown in Fig. 2. Simulations for July 1979 and for the year 1979 
represent tests of the model for two different time scales. Gaussian 
dispersion calculations are included in the shorter time period but not 
in the longer one. Those tests have been based on a relatively small 
and imperfect set of meteorological data. 

4.1.1 Data 

We now have a considerably larger meteorological data set for 
the year 1978. This is being processed to put it into suitable form for 
basic input data. Simulation for the year 1978 will then be carried 
out. This will include the calculation of averages on the annual time 
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trajectories is calculated using the total emissions for the group. 

A second method is to compute the geographic mean of the 
group, calculate trajectories for this location, then translate these 
trajectories to each emission point within the group for mass balance 
computations. For each method, care must be taken in grouping the 
sources. They must be close enough to each other (within one grid 
interval, 127 km) and in sufficiently similar terrain for the 
calculations to be valid. 

With the purpose of employing the second method, a tentative 
grouping of the 63 sources into 21 clusters was set up. However, 
based on a preliminary test of trajectory translation, it was decided 
not to group the source points for the year 1979, but to compute 
trajectories for all sources separately. Following this, tests of 
sensitivity to grouping can be performed using all the 1979 
trajectories. 



precipitation value which determines k w and k w is the next one in 
time, based on the same reasoning. 

3.2.2 Emission Inventory 

The emission inventory is the one used by the statistical model 
(Venkatram et al; 1981). The location and source strengths of 83 
sources are shown in Fig. 1. These are annual emissions. The data is 
a composite of information from several sources. U.S. emission data 
was compiled from EPA point source inventories (Benkowitz; 1979) 
and SURE data from the GCA (consulting company contracted by 
EPR1) major point and area source records. Ontario points were 
taken from the Ontario Ministry of the Environment sulphur 
emissions inventory. For the rest of Canada, emission points came 
from the AES preliminary point and area source inventory (Voldner et 
al; 1980). Area sources are treated as point sources representing the 
area. 

Trajectory and mass balance calculations are made by the 
model for 62 of the 83 sources. The remaining 21 sources lie outside 
the grid array or are too near the boundary for accurate trajectory 
calculations. Errors in the polynomial calculations occur near the 
boundaries because of the lack of data for polynomial fitting. 
Concentration and depositions resulting from the 21 sources as 
calculated by the statistical model are used as background values for 
the one year simulation. 

Some authors report grouping of sources and calculating 
representative trajectories for each group as a means of economizing 
on computer costs. One method (MOE statistical model) is to 
compute an emissions-weighted mean location for each group, then 
calculate trajectories from this location. The mass balance along the 
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less than 2% per hour while the pollutant is becoming mixed in the 
boundary layer. Smith (1981) suggests that during the initial period 
the amount of SO 4= in the puff may be in the range 5-30% of the 
total sulphur. The high value of ^used here reflects only a monthly 
value. The long term average value may be less than this. For the 
one year simulation/^ =^ % * The value ofo£ is zero. 

3.2 Data and Calculations 

Input data for the mass balance calculations consist of 
trajectories as computed in section 2, annual emissions for all source 
points and precipitation data. 

3.2.1. Precipitation Data and Scavenging Calculations 

The 1979 precipitation data which we are using was obtained 
from AES, objectively analysed to the 127 km grid. Precipitation 
amounts are 2k hour accumulations reported every 12 hours. Where 
data at a grid point was reported missing we have interpolated in 
time if (a) there was data at twelve hours before and after the 
missing time and (b) the two precipitation values were both non-zero 
or zero. This means that no interpolation was done if rain occurred 
before and no rain occurred after the missing time, or vice versa. 

For a segment end point which occurs at a precipitation 
reporting time, the value of the precipitation given at that time 
determines the value of k^ and k# for the trajectory segment. This 
approach is taken because of the nature of precipitation data. It is 
the accumulation over the preceding 24 hours. The mass of pollutant 
has been travelling through this grid square for all or part of the 
preceding 6 hours. This 6 hour interval occurs within the 2k hour 
period for which the precipitation is reported. For a segment end 
point which occurs between precipitation reporting times, the 






dissolved SO2, have been given by Scott (1981) for winter conditons. 
The coefficients are widely variable, but some order emerges when 
he classifies the precipitation as rimed or unrimed snow crystals. 
Rimed crystals grow primarily by accretion of cloud water droplets 
on the snow crystal. Unrimed crystals grow by the deposition of 
water vapour on the snow crystal, when there is /ery little liquid 
water in the cloud. The washout rate for rimed (wet) crystals is 
higher than for unrimed (dry) crystals by a factor of 50 or more. This 
compares favorably with Summers (1974) values for wet and dry snow 
which differ by a factor of 60. Scott's work was based on data from 
15 winter storms at Muskegon Michigan during the months from 
November to March. He deduced from the observations that about 
one half of these consisted of rimed snow and one half, unrimed snow. 
This result gives some support to our assumption of equal occurrence 
of wet and dry snow. 

3.1.4 Chemical Conversion Rate 

The linear rate of conversion of SO2 to SOtr in the model is an 
oversimplification of a complex process which varies with humidity, 
the presence of other contaminants and a photochemical effect. 
Different values of k t were used for the one month and the full year 
calculations. A long term average for the year was assumed to be 1% 
per hour (Newman, 1980). For a one month simulation in summer, 
when the light intensity is high, k t =2% per hour. 

Similarly, the initial conversion of SO2 to SO4, is given 
different values for the July simulation and for the full year. The 
value ^3= 10% takes into account the rapid oxidation due to high 
concentrations of SO2 and other chemicals during emission, plus 
increased photochemical effects in summer. For the initial time step 
of 6 hours this represents an increase in oxidation rate of somewhat 
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3.1.3 Wet Deposition 

Wet deposition factors k^ and k^ are loss rates which vary with 
rain intensity I. For SO2, the value k^ = 3. x 10-5 I was taken from 
Summers (197*0. This includes both below -cloud washout and in-cloud 
conversion of SO2 to SO^ = with subsequent rain-out. The in-cloud 
conversion is limited by the rate at which SO2 is drawn into the 
clouds, a value estimated (Venkatram et al;1981) as of the order 10 - ^ 
s - l. Washout of SO^= in rain as a long term value is generally 
considered to be an order of magnitude higher than for S02« The 
value \<w = 3. x 10 - ^ I was used. Further discussion of these 
relationships is given in section 5.2 

3.1.3.1 Deposition in Rain and Snow 

Summers* value of k^ for dry snow is an order of magnitude 
lower than for rain. For wet snow, however, k^ is about 6 times the 
value for rain, hence a factor of 60 times the value for dry snow. In 
all cases the snowfall rate is expressed as equivalent rainfall rate. In 
the model, for preliminary testing of a full year's data, a single value 
of kw represents conditions for the year for all grid locations. Snow 
can occur from November through March over most of the 
geographical area covered by the array. In the Great Lakes basin, 
wet snow occurs frequently. Cold dry air flowing over the relatively 
warmer open water of the lakes in winter gains heat and moisture in 
the lowest 2 to 3 km. Thus in January and February, while dry snow 
conditions may prevail north and west of the lakes, periods of wet 
snow often occur to the south and east of the individual lakes. In 
simulating a year's data, wet and dry snow conditions were assumed 
to occur with about equal frequency in winter. The average value of 
k^ for wet and dry snow is approximately equal to k# for rain. 
Hence Summers' rain values for k^ were adopted for the annual 
simulation. 

Measured values of washout coefficients for sulphate, including 



The reason for the different averaging methods for dry and wet 
depositon is that k w and k w vary spatially and temporally, while k(j 
and k(j are each constant and single-valued. 

3.1 Deposition Factors 

Empirical values from the literature have been assigned to the 
k factors of deposition and transformation. These are given in Table 
1. Each of the k factors represents a simplification of complex 
physical and chemical processes. Although these are not completely 
understood, their effects over a long period have been evaluated 
empirically. 

3.1.2 Dry Depositon 

Dry deposition is calculated as 

krf = Vj and k<j = V^ 



H H 

where V^ and V^ are deposition velocities of SO2 and SOzp 
respectively. Dry deposition depends on the type of surface and the 
atmospheric stability. The value V,j=l cm s _i was chosen for the 
model. Sehmel's (1980) literature survey shows measured values of 
V(j, the dry deposition velocity of SO2, under varying conditions of 
atmospheric stability over grass, water, snow and other surfaces as 
ranging between extremes of 0A to 7.5 cm s _1 , with the majority of 
values lying within the range 1.0 + 0.5 cm s" 1 . Chamberlain's (1980) 
overview of measured values showed the same wide range as Sehmel, 
(0.1 to 8 cm s - l) but a somewhat lower representative value, 0.4 to 
0.8 cm s" 1 . Dry depositon of SO^= is less than that of SO2 (e.g. 
Fisher; 1978). The value V c |=0.1 cm s" 1 is used in the model. The 
mixed layer depth H is constant at 1000 m, representing a long term 
average value. Possible modifications of the parameters is discussed 
in section 5.1. 



individually and not over the whole trajectory. Integrating over the 
interval t to (t + A t), 

q(t + A t) = q(t)exp(-Ki At) (3) 

s(t + A t) = s(t) exp (-K 2 A t) + q (t).3/2 k t exp (- K 2 A t). 

2(Ki-K 2 ) " 

(l-exp(-(K 1 -K 2 )At) W 

At t=0, q(0) = £ ( 1-ol-B ) (5) 

H 

and s(0) = 3/2 $ Q (6) 

H 

In eqs. (3) and (4) 

K l = k d + k w + k t 
and K 2 = k d + k w 

Q is the emission of sulphur per unit area, H is the mixed layer 
height. The factor d is the fraction of sulphur dioxide deposited 
locally (in the same grid square as the emission point) and $ 
is the fraction of emission occurring as sulphate. Equations (3) and 
(4) are computed at each segment end point of all the trajectories. 

Concentration and deposition are calculated for any particular 
area (grid square) by assuming that each puff acts independently 
from all other puffs both from the same source and from other 
sources. Thus the concentration in a grid square is a simple 
arithmetic mean of the concentrations for all trajectory segment end 
points occurring in the grid square. The average is calculated over 
the total period of time of the trajectory calculations. Such 
simplified calculations are possible because of the linear relationships 
in the equations. The average dry deposition for each grid square is 
calculated as the product of the dry deposition parameter and the 
average concentration. Wet depositon is calculated for each grid 
square every time interval then averaged over the total time period. 



-* 



3. Mass Balance 

Continuous emission from each source is modelled as a series of 
discrete puffs. Each puff is assumed to fill instantaneously a grid 
volume 127 x 127 km^ horizontally from the surface to the top of the 
mixing layer. Each puff then travels along a trajectory as calculated 
in section 2. 

Concentrations of sulphur dioxide, q and sulphate, s along the 

trajectory are calculated by linear equations which express 

concentration as a decay function. 

For sulphur dioxide, 

d^ = - (k d + k w + k t ) q (1) 

dt 



For sulphate 



ds = - (k(j + k w ) s + 3/2 k t q. 
dt 



Here k<j and k<j are dry deposition factors for sulphur dioxide and 
sulphate respectively. Similarly k w and k w are wet deposition 
factors. Chemical transformation of sulphur dioxide to sulphate is 
calculated using the factor kf This is expressed as a loss in equation 
(1) and a gain in equation (2). The proportion 3/2 is the ratio of the 
molecular weight of the two sulphur forms. This type of mass 
balance calculation is used in many models e.g. Eliassen and 
Saltbones (1975) and Bhumralkar et al (19&$. 

For a simple treatment, the k parameters are assumed 
independent of the concentration. Hence equations (1) and (2) are 
linear and may be integrated analytically. Since the value of k w and 
k w can change from one time step to another depending on whether 
rain is occurring, the analytical solution applies at each time step 



speed is reduced by 10% and the backing angle is 20°; for unstable 
conditions, the wind speed is geostrophic and the turning angle is 10°. 
Unstable conditions are assumed to prevail in the daytime and stable 
conditions at night. Discussion of more detailed calculations possible 
in the future is given in section 5.1. 

Each trajectory is modelled in segments by wind vectors 
calculated from the pressure field at each time step^,t. Constant 
acceleration along each segment is maintained using the predictor - 
corrector method of Pettersen (1956). This is a relaxation 
adjustment between the wind vectors at each end of the segment at a 
given time, t. A new trajectory is started at the source point at each 
time interval. For example, for six-hour intervals, each forty-eight 
hour trajectory consists of eight six- hour segments and a new 
trajectory is started from the source point every six hours. These 
calculations are carried out for every source point. 

Because of the random variation of the pressure field both 
spatially and temporally, each trajectory is unique. In this way the 
model calculates a collection of trajectories which can be regarded 
as a statistical ensemble of particle paths in a field of turbulence 
providing the set is large enough, i.e. if calculated for an extended 
period of time, at least one or two years. 

Many simplifying assumptions under ly the model, neglecting 
small scale complex variations. Consequently trajectories should not 
be considered to represent reality individually, but in total, as a long 
term average. Mass balance calculations, as described in the next 
section, are performed along each trajectory. The result is a long 
term average of pollutant distribution due to the trajectories. 
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2. Trajectory Calculations 

Horizontal trajectories are calculated within the mixed layer 
where the simulated wind speed and direction represent vertically 
averaged values within the layer. The advecting winds are 
geostrophic values computed from the surface pressure field and 
modified as to speed and direction to take into account the effect of 
surface friction and thermal stability. The reason for using surface 
pressures in preference to 850 mb values for the wind calculations as 
is often done, is that the time and space scales of pressure 
observations at the surface are the finest scale available. Surface 
pressures are reported hourly and six hourly at stations located 
roughly 100 km to 200 km apart. In contrast, 850 mb pressure 
observations are made every twelve hours at stations located 300 to 
1000 km apart. 

In order to create a pressure field which varies smoothly in 
space, orthogonal polynomials are fitted to the pressure data using a 
least squares fit criterion as given by Sykes and Hatton (1976). 

In the boundary layer, surface friction is of comparable 
magnitude with pressure gradient and Coriolis forces. Hence 
boundary layer winds maintain a balance among these three forces. 
The geostrophic relation, while convenient for computing winds, 
takes into account the balance between only two forces, the pressure 
gradient and Coriolis. Surface friction reduces the wind speed from 
the geostrophic value. This causes turning in the anticlockwise 
direction from the geostrophic flow, or backing, such that the wind 
flows down the pressure gradient. The magnitude of these 
modifications due to surface friction depends on the turbulence in 
the boundary layer. In the model, for stable conditions, the wind 



Thus the Lagrangian model can give results for shorter time periods 
such as a season, month, or possibly episode (4-5 days). The 
Lagrangian model traces the life history of each 6 hour 'packet' of 
emission from every source through the interactions along the flow, 
to its deposition on the earth's surface. However, the model is 
subject to uncertainties in wind computation and gross simplification 
of the physical and chemical influences during transport. Hence the 
effects must be computed as an average over many trajectories. 
Nevertheless, the averaging period can be as short as a few days if a 
sufficient number of trajectories are used and if suitable paramatric 
values are known. 

The Lagrangian model, being much more complex than the 
Statistical model, is more expensive. It takes longer to develop, 
requires a large meteorological data set and more computer time for 
its operation. In return, the Lagrangian model can compute monthly 
or episode values which the Statistical model cannot do. It can 
compute year to year variations for each month, thereby showing 
trends in concentration and deposition values at different locations. 
The Lagrangian model can also build up a bank of trajectory data to 
provide dispersion statistics needed by the Statistical model. 

In the report which follows, trajectory calculations are 
described in section 2, while mass balance of sulphur is detailed in 
section 3. The current state of the model is given in section 4, with 
projections of its operation for the near future. Results of 
preliminary tests are discussed, including a one month simulation for 
July 1979 and a simulation for the full year 1979. Recommendations 
for future modification to the model are made in section 5, with 
discussion of boundary layer winds, wet and dry depositions and 
reference to vertical dispersion and boundary layer height. 



since sulphur ranks first in concern and has parameters which are 
probably the best known and most readily computed. 

For long range transport, detailed fields near the source are not 
simulated. Each puff of pollutant is considered to be dispersed 
instantaneously in the vertical through the mixing layer. For 
calculation of a long term average (1-2 years) horizontal dispersion is 
simulated by the spread of the set of trajectories viewed as a whole 
rather than along individual trajectories. However, for simulation on 
a shorter time scale, such as a month, an episode or a day, there are 
not enough trajectories to create this dispersion. Hence it is 
necessary to apply horizontal disperson to the model calculations for 
simulation of short time periods. 

The model is applied over a section of eastern North American 
centred on the Great Lakes, on a grid array of 23 x 26 points with a 
horizontal grid distance of 127 km. Currently the duration of 
trajectories is forty-eight hours. 

The theory and development of this model has already been 
described (Atmospheric Model Development Unit; 1980), with a 
preliminary application to a sulphate episode period. 
1.1 Comparison with Statistical Model 

The Statistical model (Venkatram et al, 1981) is a simple but 
efficient model which calculates long term averages based on 
statistics of wind flow, chemical transformation and deposition. To 
compute annual average ambient concentration and depositions of 
pollutants, the model assumes a single average transport wind with a 
typical dispersion pattern. Statistical values of dry and wet duration 
times determine deposition during transport. In contrast, the 
Lagrangian model uses observed atmospheric pressures and 
precipitation at each time interval (e.g. every 6 hours) for many 
locations within the area of concern, i.e., eastern North America. 



1. Introduction 

The transport of pollutants over long distances through the 
atmosphere has aroused concern because the depositon of pollutants 
and the ambient air quality far downstream may reach harmful levels 
when the effects of many sources are added together. 

In this model of long range transport, the distances for which 
simulation is calculated are of the order 100 to 1000 km horizontally 
and one km vertically. Continuous emissions from a point source are 
modelled as a series of discrete puffs travelling successively along a 
wind flow trajectory which changes in time and space during the 
travel period. Chemical transformations of the pollutant and mass 
decay by loss to the earth's surface under wet and dry conditions are 
calculated along each trajectory. Summation of these processes over 
all trajectories for all emission sources yields concentration and 
deposition statistics at locations downstream. 

The trajectory simulations constitute the core of the model. 
Over an extended period of time, one or two years, these trajectories 
may be considered to form a statistical ensemble representing long 
term average flow patterns for the pollutant. 

The general form of the model permits the computation of mass 
budgets for any pollutant if the chemistry and the loss parameters 
can be ascertained. Currently the sulphur balance is being simulated 
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